1. Описание признаков присутствует в сопроводительном файле. Прочитаем данные. Среди данных не может быть отрицательных значений, NA обозначаются как -999.
library(readxl)
library(dplyr)
library(tidyr)
df <- read_excel("Sleep/SLEEP_shortname.xls")
df[df < 0] <- NA
head(df)
## # A tibble: 6 × 11
## NAME BODY_…¹ BRAIN…² SLOWW…³ PARADOX SLEEP LIFES…⁴ GESTT…⁵ PRED_…⁶ EXP_IND
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 African… 6654 5712 NA NA 3.3 38.6 645 3 5
## 2 African… 1 6.6 6.3 2 8.3 4.5 42 3 1
## 3 Arctic … 3.38 44.5 NA NA 12.5 14 60 1 1
## 4 Arc. gr… 0.92 5.7 NA NA 16.5 NA 25 5 2
## 5 Asian e… 2547 4603 2.1 1.8 3.9 69 624 3 5
## 6 Baboon 10.6 180. 9.1 0.7 9.8 27 180 4 4
## # … with 1 more variable: DANG_IND <dbl>, and abbreviated variable names
## # ¹BODY_WEI, ²BRAIN_WE, ³SLOWWAVE, ⁴LIFESPAN, ⁵GESTTIME, ⁶PRED_IND
## # ℹ Use `colnames()` to see all variable names
Количество строк в таблице.
nrow(df)
## [1] 62
2. Признаков немного, поэтому рассматривать будем все.
3. Кроме индексов-факторов, все признаки количественные. Индексы - порядковые признаки. Все количественные признаки непрерывные, но можно заметить дискретизацию при округлении. Проверим это, посмотрев на частоты мод и повторов.
mode_rate <- function(x) {
x <- x[!is.na(x)]
u <- unique(x)
tab <- tabulate(match(x, u))
max(tab) / length(x)
}
repeat_rate <- function(x) {
x <- x[!is.na(x)]
u <- unique(x)
(length(x) - length(u)) / length(x)
}
Отношение частоты моды к числу элементов:
df %>% summarise(BODY_WEI = round(mode_rate(BODY_WEI), 3), BRAIN_WE = round(mode_rate(BRAIN_WE), 3),
SLOWWAVE = round(mode_rate(SLOWWAVE), 3), PARADOX = round(mode_rate(PARADOX), 3),
SLEEP = round(mode_rate(SLEEP), 3), LIFESPAN = round(mode_rate(LIFESPAN), 3),
GESTTIME = round(mode_rate(GESTTIME), 3))
## # A tibble: 1 × 7
## BODY_WEI BRAIN_WE SLOWWAVE PARADOX SLEEP LIFESPAN GESTTIME
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.032 0.032 0.062 0.06 0.052 0.052 0.052
Отношение частоты повторных элементов к числу элементов:
df %>% summarise(BODY_WEI = round(repeat_rate(BODY_WEI), 3), BRAIN_WE = round(repeat_rate(BRAIN_WE), 3),
SLOWWAVE = round(repeat_rate(SLOWWAVE), 3), PARADOX = round(repeat_rate(PARADOX), 3),
SLEEP = round(repeat_rate(SLEEP), 3), LIFESPAN = round(repeat_rate(LIFESPAN), 3),
GESTTIME = round(repeat_rate(GESTTIME), 3))
## # A tibble: 1 × 7
## BODY_WEI BRAIN_WE SLOWWAVE PARADOX SLEEP LIFESPAN GESTTIME
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.032 0.048 0.188 0.4 0.241 0.19 0.155
Будем считать данные дискретными если отношение частоты моды к числу элементов больше 10%, а отношение частоты повторных элементов к числу элементов больше 15%. Тогда BODY_WEI, BRAIN_WE будут непрерывными, а BRAIN_WE, PARADOX, SLEEP, LIFESPAN, GESTTIME дискретными.
4. Не актуально для текущих данных.
5. Посмотрим на данные.
library(ggplot2)
library(GGally)
df %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"),
columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
6. Преобразуем данные. Переведем вес животного в граммы, продолжительность жизни в дни. Логарифмируем данные, чтобы было легче наблюдать корреляции (выше заметны сильно отличающиеся индивиды, это слоны). Факторизуем индексы.
dfNew <- df %>% mutate(PRED_IND = as.factor(PRED_IND), EXP_IND = as.factor(EXP_IND),
DANG_IND = as.factor(DANG_IND),
BODY_WEI = log(BODY_WEI * 1000), BRAIN_WE = log(BRAIN_WE),
LIFESPAN = log(LIFESPAN * 365.25), GESTTIME = log(GESTTIME))
Посмотрим на новые данные. Сгруппируем индивидов по факторам.
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"),
columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=PRED_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=EXP_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=DANG_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
7. Аутлайнеров нет.
8. Неоднородности нет.
9. Не актуально для текущих данных.
10. Используем описательные статистики для распределений признаков в новой выборке.
library(moments)
characteristics <- function(x) {
c(mean = round(mean(x, na.rm = TRUE), 3),
median = round(median(x, na.rm = TRUE), 3),
var = round(var(x, na.rm = TRUE), 3),
skewness = round(skewness(x, na.rm = TRUE), 3),
kurtosis = round(kurtosis(x, na.rm = TRUE) - 3, 3))
}
data.frame(lapply(as.list(dfNew[2:8]), characteristics)) %>% t()
## mean median var skewness kurtosis
## BODY_WEI 8.245 8.114 9.754 0.149 -0.441
## BRAIN_WE 3.140 2.848 5.985 0.045 -0.540
## SLOWWAVE 8.673 8.350 13.443 0.289 -0.334
## PARADOX 1.972 1.800 2.081 1.409 1.976
## SLEEP 10.533 10.450 21.222 0.196 -0.567
## LIFESPAN 8.494 8.613 0.890 -0.167 -0.841
## GESTTIME 4.444 4.360 1.124 0.049 -1.075
1. Проверим распределения признаков на нормальность.
library(reshape)
Изобразим гистограммы признаков на фоне плотности нормального распределения.
dfMelt <- dfNew %>%
dplyr::select(-PRED_IND, -DANG_IND) %>%
as.data.frame(dfNew) %>%
melt(id.vars = c("NAME", "EXP_IND"))
ggplot(dfMelt, aes(x = value)) +
geom_histogram(aes(y = ..density..)) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
geom_line(aes(y = dnorm(value,
mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
)), color = "blue")
Изобразим выборочную функцию распредления признаков на фоне функции распределения нормального распределения.
ggplot(dfMelt, aes(x = value)) +
stat_ecdf() +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
geom_line(aes(y = pnorm(value,
mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
)), color = "blue")
Давайте еще нарисуем PP-plot и QQ-plot.
library(qqplotr)
ggplot(dfMelt, aes(sample = value)) +
stat_pp_point(size = 1) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_pp_line(color = "blue")
ggplot(dfMelt, aes(sample = value)) +
stat_qq_point(size = 1) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_qq_line(color = "blue")
2. Похоже, что нормальность имеет место. Проверим это с помощью критериев.
library(nortest)
У нас достаточно повторений в некоторых выборках, чтобы не использовать критерий Пирсона, но все ранво проведем его, скорее, ради академического интереса. Критерий Лиллиефорса - это модификация критерия Колмогорова-Смирнова для проверки сложных гипотез нормальности данных. Критерий Андерона-Дарлинга - это один из критериев типа w^2. Критерий Шапира-Уилка - примерно квадрат корреляции между x и y в n.p.p.
normality.tests <- function(x) {
c(Pearson = round(pearson.test(x)$p, 5),
Lilliefors = round(lillie.test(x)$p, 5),
Anderson.Darling = round(ad.test(x)$p, 5),
Shapiro.Wilk = round(shapiro.test(x)$p, 5))
}
data.frame(lapply(as.list(dfNew[2:8]), normality.tests)) %>% t()
## Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## BODY_WEI 0.26277 0.04414 0.46210 0.72717
## BRAIN_WE 0.71196 0.70193 0.64143 0.69703
## SLOWWAVE 0.86095 0.78346 0.90556 0.77036
## PARADOX 0.13013 0.02883 0.00027 0.00009
## SLEEP 0.07080 0.81866 0.37487 0.13963
## LIFESPAN 0.47137 0.45797 0.20207 0.29791
## GESTTIME 0.51083 0.22157 0.11163 0.12164
Гипотеза о нормальности распределения для быстрого сна отвергается многими критериями с уровнем значисости 0.05. Гипотеза нормльности распределения для веса животного отвергается критерием Лиллиефорса с уровнем значисости 0.05. Остальные критерии гипотезы не отвергли с уровнем значисости 0.05.
3. Опишем разницу между животными по степени опасности места, где они спят.
dfNew %>% ggplot(aes(y=BODY_WEI, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=BRAIN_WE, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=SLOWWAVE, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=PARADOX, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=SLEEP, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=LIFESPAN, colour=EXP_IND)) + geom_boxplot()
dfNew %>% ggplot(aes(y=GESTTIME, colour=EXP_IND)) + geom_boxplot()
Можно заметить, что животные, живущие в более защищенных местах, 1,2) имеют меньший вес тела и мозга; 3,4,5) чаще дольше спят в общем и в разных фазах сна по отдельности; 6,7) имеют менее длительные продолжительность жизни и период вынашивания потомства.
Давайте проверим наблюдения 1,3) c помощью критериев для животных, живущих в ниболее и в наименее защищенных местах.
dfNew1 <- dplyr::filter(dfNew, EXP_IND == 1)
dfNew5 <- dplyr::filter(dfNew, EXP_IND == 5)
nrow(dfNew1)
## [1] 27
nrow(dfNew5)
## [1] 13
Ранее мы узнали, что многие признаки распределены нормально, но сейчас мы будем сравнивать разные группы индивидов, признаки внутри которых могут уже не имеют нормальность.
ggplot(dfNew1, aes(sample = BODY_WEI)) +
stat_qq_point(size = 1) +
labs(x = "", y = "") +
stat_qq_line(color = "blue")
ggplot(dfNew1, aes(sample = SLEEP)) +
stat_qq_point(size = 1) +
labs(x = "", y = "") +
stat_qq_line(color = "blue")
ggplot(dfNew5, aes(sample = BODY_WEI)) +
stat_qq_point(size = 1) +
labs(x = "", y = "") +
stat_qq_line(color = "blue")
ggplot(dfNew5, aes(sample = SLEEP)) +
stat_qq_point(size = 1) +
labs(x = "", y = "") +
stat_qq_line(color = "blue")
Используем те же критерии, что использовали ранее.
Вес тела для животных, спящих в более защищенном месте.
normality.tests(dfNew1$BODY_WEI)
## Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## 0.07419 0.42187 0.30137 0.24588
Сон для животных, спящих в более защищенном месте.
normality.tests(dfNew1$SLEEP)
## Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## 0.16492 0.16347 0.09438 0.12591
Вес тела для животных, спящих в менее защищенном месте.
normality.tests(dfNew5$BODY_WEI)
## Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## 0.36851 0.78725 0.91119 0.99223
Сон для животных, спящих в менее защищенном месте.
normality.tests(dfNew5$SLEEP)
## Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## 0.00889 0.00047 0.00476 0.00526
Нормальность отверглась с уровням значимости 0.05 толко для сная животных, спящих в менее защищенном месте. Причем по всем критериям отверглась.
Итак. Группы у нас независимые и некоторые из них имеют нормальное распределение.
4. Сравним распределения с помощью t-критериев для независимых выборок. Но сначала проверим гипотезу о равенстве дисперсий распределений.
Для веса тела можно использовать критерий Фишера, а для сна нельзя, так как для сна не выполняется нормальность, поэтому ему доверять не будем, но проведем ради академического интереса. Так же используем критерий Левена с модулями.
var.tests <- function(x, y) {
c(Fisher = round(var.test(x, y)$p.value, 5),
Levene = round(t.test(abs(x - mean(x, na.rm = TRUE)),
abs(y - mean(y, na.rm = TRUE)))$p.value, 5))
}
Вес тела:
var.tests(dfNew1$BODY_WEI, dfNew5$BODY_WEI)
## Fisher Levene
## 0.44453 0.28662
Сон:
var.tests(dfNew1$SLEEP, dfNew5$SLEEP)
## Fisher Levene
## 0.01559 0.00497
Гипотеза о равных дисперсиях не отверглась для веса, а для сна отверглась с уровнем значимости 0.05.
Посмотрим на t-критерии с равными и неравными дисперсиями для независимых выборок.
t.tests <- function(x, y) {
c(Two.sample.t.test = t.test(x, y, var.equal=TRUE)$p.value,
Welch.two.sample.t.test = t.test(x, y, var.equal=FALSE)$p.value)
}
Вес тела:
t.tests(dfNew1$BODY_WEI, dfNew5$BODY_WEI)
## Two.sample.t.test Welch.two.sample.t.test
## 6.044937e-07 4.735032e-07
Сон:
t.tests(dfNew1$SLEEP, dfNew5$SLEEP)
## Two.sample.t.test Welch.two.sample.t.test
## 1.292774e-07 1.676559e-10
Первый критерий точный, второй используется не по назначению (так как дисперсии разные), третий и четвертный ассимптотические. Как и предполагалось отверглись гипотезы в пользу того, что мы наблюдали на box-plot-ах выше.
5,6. Так как в случае сна удалось проверить гипотезу только асимптотическим критерием, а выборки у нас небольшие, воспользуемся другими критериями. Выбросов у нас нет, но все же используем критерий Уолкоксона ради академического интереса.
Вес:
wilcox.test(dfNew1$BODY_WEI, dfNew5$BODY_WEI, paired=FALSE)$p.value
## [1] 2.183504e-05
Сон:
wilcox.test(dfNew1$SLEEP, dfNew5$SLEEP, paired=FALSE)$p.value
## [1] 1.002996e-05
Мощность у этого критерия меньшая, но гипотезы так же отверглись.
7. Наконец используем двухвыборочный критерий Колмогорова-Смирнова, который имеет меньшую мощность, но у которого более общая альтернативная гипотеза.
Вес:
ks.test(dfNew1$BODY_WEI, dfNew5$BODY_WEI)$p.value
## [1] 0.0001415934
Сон:
ks.test(dfNew1$SLEEP, dfNew5$SLEEP)$p.value
## [1] 4.413664e-05
Гипотезы так же отверглись.
8. Не актуально для текущих данных.
1. Посмотрим еще раз на ggpairs plot.
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"),
columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=PRED_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=EXP_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
dfNew %>% dplyr::select(-NAME) %>%
ggpairs(diag=list(continuous = "barDiag"), aes(colour=DANG_IND), legend=1,
columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))
Можно заметить 1) положительную корреляцию веса тела, веса мозга, продолжительности жизни и периода вынашивания потомства; 2) положительную корреляцию фаз сна и сна в общем смысле; 3) отрицательную корреляцию всего, что было перечислено в пункте 1), со всем, что было перечислено в пункте 2).
2. Посмотрим на корреляцию Пирсона между признаками.
dplyr::select(dfNew, -NAME) %>%
mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
cor(method = "pearson", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 3))) +
scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) +
ggtitle("dfNew pearson") +
theme(axis.text.x = element_text(angle = 50, hjust = 1))
Видим подтверждение того, что наблюдали на ggpairs, еще заметны отрицательные корреляции индексов и сна. То есть чем в большей опасности находится животное, тем меньше оно спит. Еще заметна сильная корреляция между размерами животного и незащищенностью места, где оно спит.
3. Теперь посмотрим на корреляции Спирмана между признаками.
dplyr::select(dfNew, -NAME) %>%
mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
cor(method = "spearman", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 3))) +
scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) +
ggtitle("dfNew spearman") +
theme(axis.text.x = element_text(angle = 50, hjust = 1))
Больших изменений нет (так как у нас нет выбросов и распределения близкие к нормальным).
4. Также посчитаем интересные частные корреляции.
library(ppcor)
Посмотрим на частную корреляцию веса животного и сна за вычетом идексов опасности и на частную корреляцию степени защищенности животного во время сна и сна за вычетом веса и связанных с ним критериев.
Это будет иметь смысл, так как между индексами и весом/сном коэффициент Пирсона значим, то есть есть линейная зависимость между индексами и весом/сном, которая будет вычитаться при подсчете частных корреляций.
Частая корреляция веса животного и сна за вычетом идексов опасности.
((dplyr::select(dfNew, -NAME, -BRAIN_WE, -LIFESPAN, -GESTTIME, -SLOWWAVE, -PARADOX) %>%
mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
drop_na() %>%
pcor(method = "spearman"))$estimate %>%
as.data.frame())["SLEEP", "BODY_WEI"]
## [1] -0.2668484
Отрицательная корреляции между сном и весом животного остается. Животные с большим весом спят меньше.
Часнтая корреляция степени защищенности животного во время сна и сна за вычетом веса и связанных с ним критериев.
((dplyr::select(dfNew, -NAME, -PRED_IND, -DANG_IND, -SLOWWAVE, -PARADOX) %>%
mutate(EXP_IND = as.numeric(EXP_IND)) %>%
drop_na() %>%
pcor(method = "spearman"))$estimate %>%
as.data.frame())["SLEEP", "EXP_IND"]
## [1] -0.4055336
Отрицательная корреляция между сном и индексом защищенности животного во время сна остается. Животные, находящиеся в большей опасности во время сна, спят меньше.